Post-collapse dynamics of self-gravitating 
Brownian particles and bacterial populations 

Clement Sire and Pierre-Henri Chavanis 
February 2, 2008 



Laboratoire de Physique Theorique (FRE 2603 du CNRS), Universite Paul Sabatier, 
118, route de Narbonne, 31062 Toulouse Cedex 4, France 
E-mail: Clement.Sire@irsamc.ups-tlse.fr & Chavanis@irsamc.ups-tlse.fr 



Abstract 

We address the post-collapse dynamics of a self-gravitating gas of Brownian parti- 
cles in D dimensions, in both canonical and microcanonical ensembles. In the canonical 
ensemble, the post-collapse evolution is marked by the formation of a Dirac peak with in- 
creasing mass. The density profile outside the peak evolves self-similarly with decreasing 
central density and increasing core radius. In the microcanonical ensemble, the post- 
collapse regime is marked by the formation of a "binary" -like structure surrounded by an 
almost uniform halo with high temperature. These results are consistent with thermody- 
namical predictions in astrophysics. We also show that the Smoluchowski-Poisson system 
describing the collapse of self-gravitating Brownian particles in a strong friction limit is 
isomorphic to a simplified version of the Keller-Segel equations describing the chemotactic 
aggregation of bacterial populations. Therefore, our study has direct applications in this 
biological context. 

1 Introduction 

Self-gravitating systems such as globular clusters and elliptical galaxies constitute a Hamilto- 
nian system of particles in interaction that can be supposed isolated in a first approximation 
PQ. Since energy is conserved, the proper description of stellar systems is the microcanonical 
ensemble j2]. The dynamical evolution of elliptical galaxies is governed by the Vlasov-Poisson 
system which corresponds to a collisionless regime. On the other hand, the kinetic theory of 
globular clusters is based on the Landau-Poisson system (or the orbit averaged Fokker-Planck 
equation) which describes a collisional evolution. These equations conserve mass and energy. 
Furthermore, the Landau equation increases the Boltzmann entropy (H-theorem) due to stellar 
encounters. These equations have been studied for a long time in the astrophysical literature 
and a relatively good physical understanding has now been achieved [T]. In particular, globu- 
lar clusters can experience core collapse j3 El E] related to the "gravothermal catastrophe" 
concept [Jj. 

For systems with long-range interactions, statistical ensembles are not equivalent j2] • There- 
fore, it is of conceptual interest to compare the microcanonical evolution of stellar systems to a 
canonical one in order to emphasize the analogies and the differences. This can be achieved by 
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considering a gas of self-gravitating Brownian particles [S] subject to a friction originated from 
the presence of an inert gas and to a stochastic force (modeling turbulent fluctuations, colli- 
sions,...). This system has a rigorous canonical structure where the temperature T measures 
the strength of the stochastic force. Thus, we can precisely check the thermodynamical predic- 
tions of Kiessling and Chavanis ^U] obtained in the canonical ensemble. In the mean-field 
approximation, the self-gravitating Brownian gas model is described by the Kramers-Poisson 
system. In a strong friction limit (or for large times) it reduces to the Smoluchowski-Poisson 
system. These equations conserve mass and decrease the Boltzmann free energy [TT]. They 
possess a rich physical and mathematical structure and can lead to a situation of "isothermal 
collapse", which is the canonical version of the gravothermal catastrophe. These equations 
have not been considered by astrophysicists because the canonical ensemble is not the correct 
description of stellar systems and usual astrophysical bodies do not experience a friction with 
a gas (except dust particles in the solar nebula [12J). Yet, it is clear that the self-gravitating 
Brownian gas model is of considerable conceptual interest in statistical mechanics to understand 
the strange thermodynamics of systems with long-range interactions and the inequivalence of 
statistical ensembles. In addition, it provides one of the first model of stochastic particles with 
long-range interactions, thereby extending the classical Einstein-Smoluchowski model ^3] to a 
more general context [TT] . 

In addition, it turns out that the same type of equations occur in biology in relation with the 
chemotactic aggregation of bacterial populations [Tlj. A general model of chemotactic aggre- 
gation has been proposed by Keller Sz Segel ^5] i n the form of two coupled partial differential 
equations. In some approximation, this model reduces to the Smoluchowski-Poisson system [8] . 
Therefore, there exists an isomorphism between self-gravitating Brownian particles and bacte- 
rial colonies. Non-local drift- diffusion equations analogous to the Smoluchowski-Poisson system 
have also been introduced in two-dimensional hydrodynamics in relation with the formation of 
large-scale vortices such as Jupiter's great red spot [TBJ [T71 [TH] . These analogies give further 
physical interest to our Brownian model. 

In a recent series of papers [HJ [20] , we have studied the dynamics and thermodynamics 
of self-gravitating Brownian particles confined within a spherical box of radius R in a space 
of dimension D. In these works, we focused on the pre-collapse regime. In the canonical 
situation (fixed temperature T) we showed that a critical temperature T c exists below which 
the system undergoes a gravitational collapse leading to a finite time singularity at t — t co u. 
For t — > t coiiy the evolution is self-similar in the sense that the density profile evolves as p(r, t) = 
Po(t) f (r I r (t)) where f(x) is independent on time. For x — > +oo, f(x) ~ x~ a with a = 2. The 
central density increases as po ~ (tcoii — t)^ 1 and the core radius decreases as r ~ (t co u — t) l l a . 
For T = 0, the exponent is a = 2D/(D + 2). The scaling profiles can be calculated analytically. 
These results can be compared with those of Penston [2JJ who considered an isothermal collapse 
modelled by the Euler- Jeans equations. We also introduced a microcanonical description of 
Brownian particles by letting the temperature T(t) evolve in time so as to conserve energy. 
This can provide a simplified model for the violent relaxation of collisionless stellar systems 
[TT] or, simply, a numerical algorithm ^JJ to determine what is the maximum entropy state 
at fixed E and M. In the microcanonical situation, there exists a critical energy E c (Antonov 
energy) below which the system collapses. For t — > t co u, there exists a pseudo-scaling regime 
where a passes very slowly from a max = 2.21... to a = 2. Numerical simulations suggest that 
T(t) remains finite at t = t co u so that the true scaling regime corresponds to a = 2, as in the 
canonical situation (201 122] • 

What happens after t co u ? By investigating the case T = 0, we found in [THj that the 
evolution continues in the post-collapse regime with the formation of a Dirac peak accreting 
more and more mass as M{t) ~ (t — t co u) D l 2 while the density outside the peak evolves self- 
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similarly with decreasing central density po ~ (t — tcoii) 1 an d increasing core radius r ~ 

D-\-2 

(t — tcoii)^ m ~ ■ Our aim in this paper is to investigate the post-collapse regime for T 7^ in 
both canonical and microcanonical ensembles. In Sec. we emphasize the analogy between 
self-gravitating Brownian particles and bacterial populations. In the following, we shall use 
the astrophysical terminology but we stress that our results apply equally well to biology 
where the chemotactic model has more concrete physical applications. In Sec. El we set the 
notations and recall the main results concerning the pre-collapse dynamics. In Sec. IH we 
study the post-collapse dynamics at T = by a method different from [19] . which can be 
generalized at finite temperature. The post-collapse dynamics at T > is precisely considered 
in Sec. In the canonical ensemble, we show that the system forms a Dirac peak whose mass 
increases as M(t) ~ {t — t C oii) D ^ 2 ~ l while the density profile for r > expands self-similarly with 
Po ~ (t—tcou)^ 1 and r ~ (t—tcoii) 1 ^ 2 - For large times, the system is made of a Dirac peak of mass 
~ M surrounded by a light gas of Brownian particles (with negligible self- interaction). Due to 
thermal motion, complete collapse takes an infinite time (contrary to the case T = 0). For t — > 
+00, the mass contained in the Dirac peak increases as 1 — M(t)/M ~ exp(— At) where A is the 
fundamental eigenvalue of a quantum problem. For T — > 0, we find that A = 1/4T+cd/T 1//3 -|-.... 
In the microcanonical ensemble, the post-collapse regime is very pathological. The system tends 
to create a "Dirac peak of + mass" surrounded by a uniform halo with infinite temperature. 
The central structure is reminiscent of a "binary star" containing a weak mass 2m <C M = Nm 
but a huge binding energy comparable to the potential energy of the whole cluster. Since our 
model is essentially mean-field, the physical formation of binaries is replaced by the type of 
structures mentioned above. The formation of a "Dirac peak" containing the whole mass in 
canonical ensemble and the formation of "binaries" in microcanonical ensemble are expected 
from thermodynamical considerations [231 E3 121 HOI EH El • These are the structures which 
provoke a divergence of free energy (at fixed mass and temperature) and entropy (at fixed 
mass and energy) respectively All these analytical results are confirmed by numerical 

simulations of the Smoluchowski-Poisson system in Sees. and [7| We were able in particular 
to "cross the singularity" at t — t co u and describe the post-collapse dynamics. 



2 Analogy between self-gravitating Brownian particles 
and bacterial populations 

2.1 Self-gravitating Brownian particles 

We consider a system of self-gravitating Brownian particles described by the iV coupled stochas- 
tic equations 

(1) = v 4 , = -£v, - ViUfa, r N ) + V2D' Ri(t), 

where £ is the friction coefficient, D' is the diffusion coefficient and Rj(t) is a white noise 
satisfying (Rj(t)) = and (R a>i {t)R h j{t')) = 5ij5 a bS(t — £'), where a,b = 1,...,D refer to the 
coordinates of space and i, j = 1, iV to the particles. The particles interact via the potential 
U(r\, r;v) = J2i<j u ( r i ~ r j)- I n this paper, u(rj — r 3 -) is the Newtonian binary potential in 
D dimensions. The stochastic process defines a toy model of gravitational dynamics which 
extends the classical Brownian model ^3] to the case of stochastic particles in interaction. In 
this context, the friction is due to the presence of an inert gas and the stochastic force is due 
to classical Brownian motion, turbulence, or any other stochastic effect. 
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Starting from the iV-body Fokker-Planck equation and using a mean-field approximation 
we can derive the nonlocal Kramers equation 



df df df d f ,df 
^ dt dr dv dv\ dv ^ 

where F = — V$ is the smooth gravitational force felt by the particles. The gravitational 
potential $ is related to the density p = J fd D v by the Poisson equation 

(3) A$ = S D Gp, 

where Sjj is the surface of the unit .D- dimensional sphere. Equation (J2J) can be considered as 
a generalized version of the Kramers- Chandrasekhar equation introduced in a homogeneous 
medium |27]. In this work, the diffusion and the friction model stellar encounters in a simple 
stochastic framework. The condition that the Maxwell-Boltzmann distribution is a stationary 
solution of Eq. (J2J) leads to the Einstein relation £ = D'f3 where (3 = 1/T is the inverse 
temperature (we have included the mass of the particles and the Boltzmann constant in the 
definition of T). In our case, we do not assume that the medium is homogeneous, so that we have 
to solve the Kramers-Poisson system. This makes the study much more complicated than usual. 
Up to date, we do not know any astrophysical application of this model although there could be 
connexions with the process of planetesimal formation in the solar nebula jT^j. Whatever, this 
model is interesting to develop on a conceptual point of view because it possesses a rigorous 
thermodynamical structure and presents the same features as more realistic models (isothermal 
distributions, collapse, phase transitions,...). For this Brownian model, the relevant ensemble is 
the canonical one since the temperature T is fixed. Therefore, the Kramers-Poisson system can 
be viewed as the canonical counterpart of the Landau- Poisson system. It is interesting to study 
these two models in parallel to illustrate dynamically the inequivalence of statistical ensembles 
for systems with long-range interactions. 

To simplify the problem further, we shall consider the strong friction limit £ — > +00, or 
equivalently the limit of large times t 3> In that approximation, we can neglect the inertia 
of the particles. Then, the coupled stochastic equations (JTJ) simplify to 

dri 



(4) ^ = -ViC^rx, r N ) + V2D' R,(t). 
Furthermore, to leading order, the velocity distribution is Maxwellian 

(5) /(r,v,t)= f P(r,t)e-^, 



2ir, 

and the Kramers equation reduces to the Smoluchowski equation 

"1 



< 6 > I-* 



^(TVp + pV$) 



It can be shown that the Kramers equation decreases the Boltzmann free energy 
(7) F[f] =E-TS = J /yd D rd D v + i J p$d D r + T J f\nfd D rd D w 1 

i.e. F < and F = at statistical equilibrium (canonical if-theorem) . Similarly, the Smolu- 
chowski equation decreases the free energy F[p] which is obtained from F[f] by using the fact 



4 



that the velocity distribution is Maxwellian in the strong friction limit. This leads to the 
classical expression 



(8) 



F[p] = T p\npd n r+- / p$ d v r 



The passage from the Kramers equation to the Smoluchowski equation in the strong friction 
limit is classical [2B|- It can also be obtained formally from a Chapman-Enskog expansion [29 . 

In order to prevent finite time singularities and infinite densities, we can consider a model of 
self-gravitating Brownian fermions, enforcing the constraint / < rj (Pauli exclusion principle). 
The corresponding Kramers equation takes the form 



(9) 



dt dr 



,a/ = d_ 

<9v <9v 



In the strong friction limit, we obtain a Smoluchowski equation of the form 

'1 



(10) 



dp 
di 



(Vp + pV$) 



where p(p) is the equation of state of the Fermi gas. The fermionic Smoluchowski-Poisson sys- 
tem has been studied in [3U] . Generalized Kramers and Smoluchowski equations are introduced 
in P]. 



2.2 The Keller-Segel model 

The name chemotaxis refers to the motion of organisms (amoeba) induced by chemical signals 
(acrasin). In some cases, the biological organisms secrete a substance that has an attractive 
effect on the organisms themselves. Therefore, in addition to their diffusive motion, they move 
systematically along the gradient of concentration of the chemical they secrete (chemotactic 
flux). When attraction prevails over diffusion, the chemotaxis can trigger a self-accelerating 
process until a point at which aggregation takes place. This is the case for the slime mold 
Dictyostelium Discoideum and for the bacteria Escherichia coli |14j . 

A model of slime mold aggregation has been introduced by Keller & Segel in the form 
of two coupled differential equations 

(11) % = V(D 2 Vp)-V( J D 1 Vc), 



(12) — = -k{c)c + f{c)p + D c Ac. 

at 

In these equations p(r, t) is the concentration of amoebae and c(r, t) is the concentration of 
acrasin. Acrasin is produced by the amoebae at a rate /(c). It can also be degraded at 
a rate k(c). Acrasin diffuse according to Fick's law with a diffusion coefficient D c . Amoebae 
concentration changes as a result of an oriented chemotactic motion in the direction of a positive 
gradient of acrasin and a random motion analogous to diffusion. In Eq. (jllj) . /^(p, c) is the 
diffusion coefficient of the amoebae and D\ (p, c) is a measure of the strength of the influence of 
the acrasin gradient on the flow of amoebae. This chemotactic drift is the fundamental process 
in the problem. 
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A first simplification of the Keller-Segel model is provided by the system of equations 
(13) ^ = J DAp- X V(pVc), 



dc 

(14) — = D'Ac + ap-bc, 

where the parameters are positive constants. An additional simplification, introduced by Jager 
& Lauckhaus [HI] consists in ignoring the time derivative in Eq. (|14jl. This is valid in the case 
where the diffusion coefficient D' is large. Taking also b = 0, we obtain 

(15) ^ = DAp- X V(pVc), 



(16) Ac = -Ap, 

where A = a/D'. Clearly, these equations are isomorphic to the Smoluchowski-Poisson system 
(jnj-© describing self-gravitating Brownian particles in a strong friction limit. In particular, 
the chemotactic flux plays the same role as the gravitational drift in the overdamped limit of the 
Brownian model. When chemotactic attraction prevails over diffusion, the system is unstable 
and the bacteria start to aggregate. This blow-up is similar to the collapse of self-gravitating 
systems in a canonical situation. We note that in the Keller-Segel model, the diffusion coefficient 
can depend on the density, leading to anomalous diffusion. Such a situation is considered in 
J2D] where the nonlinear Smoluchowski-Poisson system is studied. 

The Keller-Segel model ignore clumping and sticking effects. However, at the late stages 
of the blow-up, when the density of amoebae has reached high values, finite size effects and 
stickiness must clearly be taken into account. As a first step, we can propose j2H] to replace 
the classical equation ([15)1 by an equation of the form 

(17) ^ = jDA p- x V(p(o- -p)Vc), 

which enforces a limitation p(r, t) < <Jq on the maximum density of amoebae. This is the 
counterpart of the model of self-gravitating Brownian fermions [3U]- These type of non-local 
Fokker-Planck equations also occur in 2D hydrodynamics and astrophysics in relation with the 
formation of large-scale vortices and galaxies ^TJEI]. Their systematic study is clearly of broad 
interest [TTj . 



Collapse dynamics of self-gravitating Brownian parti- 
cles 



3.1 The Smoluchowski-Poisson system 

At a given temperature T controlling the diffusion coefficient, the density p(r, t) of self-gravitating 
Brownian particles satisfies the following coupled equations: 



(18) 
(19) 



dp 
dt 
A$ 



A 
SpGp 



(TVp + pV$) 
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where $ is the gravitational potential, and Sd is the surface of the unit .D-dimensional sphere. 

From now on, we set M = R = G = ^ = 1 and we restrict ourselves to spherically symmetric 
solutions. The equations of the problem become 

(20) ^ = V(TVp + pV$), 

(21) A$ = S D p, 

with proper boundary conditions in order to impose a vanishing particle flux on the surface of 
the confining sphere. These read 

(22) ^(0,t) = 0, *(1)=23d, ^(l)+p(l) = 0, 

for D > 2. For £) = 2, we take $(1) = on the boundary. Integrating Eq. f!21l) once, we can 
rewrite the Smoluchowski-Poisson system in the form of a single integrodifferential equation 

The total energy is given as the sum of the kinetic and potential contributions 

(24) E = ^T+ l -J pZ>d D r. 

The Smoluchowski-Poisson system is also equivalent to a single differential equation 

dM rT1 fd 2 M D-ldM\ 1 <9M 

(25) ir= T hj- I oT" + rr^^- 



for the quantity 

(26) M(r,t)= / p{r f )S D r' D - 1 di J t 

Jo 

which represents the mass contained within the sphere of radius r. The appropriate boundary 
conditions are 

(27) M(0,t) = N (t), Af(M) = l, 

where iVo(t) = 0, except if the density develops a condensed Dirac peak contribution at r = 0, of 
total mass N (t). It is also convenient to introduce the function s(r,t) = M(r,t)/r D satisfying 

— -T^— D+l dS \ ( ^£ 

<9t ySr 2 r dr J \ dr J 



3.2 Self-similar solutions of the Smoluchowski-Poisson system 

In |SJ El 120] > we have shown that in the canonical ensemble (fixed T) , the system undergoes 
gravitational collapse below a critical temperature T c depending on the dimension of space. 
The density develops a scaling profile, and the central density grows and diverges at a finite 
time t coU . The case D = 2 was extensively studied in [T5] and turns out to be very peculiar. 
Throughout this paper, we restrain ourselves to the more generic case D > 2, although other 
dimensions play a special role as far as static properties are concerned (see pi, 20J). 
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We look for self-similar solutions of the form 



T \ 1/2 



(29) P(r,t) = po(t)f — - , r = - 

where the King's radius r defines the size of the dense core PQ. In terms of the mass profile, 
we have 

(30) M(r 1 t) = M (t)g(-^\ with M (t) = p r D , 
and 

(31) g(x) = S D f f{x')x D - l dx'. 

Jo 

In terms of the function s, we have 

(32) s ^t) = Po (t)s(^—\ with S(x) 9 ^ 



r (t) J ' x D 
Substituting the ansatz (^2j) into Eq. (J2EI, we find that 

(33) - ^^xS'(x) = p\ ( S"(x) + ^±±S'(x) + xS(x)S'(x) + DS(x) 2 ) , 
at r at \ x J 

where we have set x = r/r . The variables of position and time separate provided that pQ 2 dp /dt 
is a constant that we arbitrarily set equal to 2. After time integration, this leads to 

(34) Po{t) = ^{t co u-t)-\ 

so that the central density becomes infinite in a finite time t co u. The scaling equation now reads 



(35) 2S + xS' = S" + ^^S' + S(xS' + DS). 

x 

The scaling solution of Eq. was obtained analytically in and reads 

4 



(36) S(x) 



D - 2 + x 2 ' 

which decays with an exponent a = 2. This leads to 

,nn\ r, ^ 4(D - 2) Z 2 + D . . 4x D 



S D (D-2 + x 2 ) 2 ' * w D-2 + a; 2 

Note finally that within the core radius tq, the total mass in fact vanishes as t — > t co «. 
Indeed, from Eq. (jSOJ), we obtain 



(38) M{r (t),t) ~ p (t)r D (t) ~ T D / 2 (t co « - t) 



D/2-1 



Therefore, the collapse does not create a Dirac peak ("black hole"). 

In [TU] . we have also studied the collapse dynamics at T = for which we obtained 



(39) po(t) ~ S^itcou-t) 



-1 
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as previously, but the core radius is not given anymore by the King's radius which vanishes for 
T = 0. Instead, we find 



(40) r ~ Po 
with 

/ \ 2D 

(41) a 



l/a 



D + 2 

The scaling function S(x) is only known implicitly 



(42) 



2 S(x) 



D + 2 



D 



D+2 _ 2R _ 

= Kx D + 2 S(x) 



where K is a known constant (see [TH] for details), 5(0) = jy^, and the large x asymptotics 
S(x) ~ f(x) ~ x~ a . The mass within the core radius is now 

(43) M(r (t),t) ~ Po (t)r D (t) ~ (t cM -t) D / 2 , 

and it again tends to zero as t — > t co u. Comparing Eq. (|S5|l and Eq. (|4*3j) suggests that if the 
temperature is very small, an apparent scaling regime corresponding to the T = case will 
hold up to a cross-over time £*, with 

(44) tcou ~ U ~ T D ^ 2 . 
Above i*, the T^O scaling ultimately prevails. 



4 Post-collapse dynamics at T = 

So far, all studies concerning the collapse dynamics of self-gravitating Brownian particles have 
concentrated on the time period t < t co u. A natural question arises: what is happening for 
t > tcou ? The first possible scenario is that the state reached at t — t co u is in fact a stationary 
state. However, its is easy to check (see jS]) that this is absolutely not the case. In addition, the 
preceding study leads to a sort of paradox |24J. Indeed, we know that the statistical equilibrium 
state in the canonical ensemble is a Dirac peak jUJ^]. This is not the structure that forms at 
t = t co u. This structure is singular at the origin (p ~ r~ 2 ) but different from a Dirac peak (in 
particular the central mass is zero). This means that the evolution must continue after t co ii- 
In particular, we will show that the Dirac peak predicted by statistical mechanics forms in the 
post-collapse regime. 

The scenario that we are now exploring is the following. A central Dirac peak containing a 
mass No(t) emerges at t > t co ii, whereas the density for r > satisfies a scaling relation of the 
form 

(45) p(r,i)=p (t)/(-^), 

where p (t) is now decreasing with time (starting from p (t = t co u) +00) and r (t) grows 
with time (starting from r (t = t co u) = 0). As time increases, the residual mass for r > is 
progressively swallowed by the dense core made of particles which have fallen on each other. It 
is the purpose of the rest of this paper to show that this scenario actually holds, as well as to 
obtain analytical and numerical results illustrating this final collapse stage. 
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In this section, we present an alternative treatment to that of where this scenario was 
analytically shown to hold at T = 0. This new approach is a good introduction to the general 
T ^ case which is studied in the next section. We refer the reader to jTH] for an explicit 
solution of the T = post-collapse regime, which we found, leads to a central peak containing 
all the mass in a finite time t en &. 

For T = 0, the dynamical equation for the integrated mass M(r,t) reads 

dM 1 dM 

with boundary conditions 

(47) M(0,t) = N (t), M(l,t) = l. 
We define po such that for small r 

(48) M(r,t) - N (t) = Po (t) r -^ + ... 

Up to the geometrical factor Sp , po(t) is the central residual density (the residual density is 
defined as the density after the central peak has been subtracted). For r = 0, Eq. flfijl leads to 
the evolution equation for iVo 

(49) f! =fl>Wo . 

As N (t) = for t < tcoii, and since this equation is a first order differential equation, it looks 
like N (t) should remain zero for t > t co u as well. However, since po(t co u) = +oo, there is 
mathematically speaking no global solution for this equation and non zero values for N (t) can 
emerge from Eq. (|49jh as will soon become clear. 
We then define 

(50) s{r,t) = , 



r 1 



which satisfies 



ds (da n \ N ( ds r , , 

By definition, we have also s(0, t) = p (t)/D. 
We now look for a scaling solution of the form 

r 



ro(t) 



(52) s(r,t) = p (t)S 
with 5(0) = D- 1 and 

(53) po(t) = r (t)- a , 

where vq is thus defined without ambiguity. Inserting this scaling ansatz in Eq. (j51|) . and 
defining the scaling variable x = t/tq, we find 



(54) (aS + xS') = S(DS + xS') + J%±-(DS + xS' - 1). 

apQ at Po r o x 
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Imposing scaling, we find that both time dependent coefficients appearing Eq. (|54|) should be 
in fact constant. We thus define a constant p such that 



(55) iV = uporj?, 
and set 

(56) — = —k, 

apQ at 

with k > 0, as the central residual density is expected to decrease. Equation (foljj) implies that 
Po ~ (t — tcou)^ 1 , which along with Eq. (jS^jl implies that N ~ (t — t co u) D l a ~ l . We thus find a 
power law behavior for Nq, which in order to be compatible with Eq. (|49|). leads to 



(57) po(t) = i^- 1 ) tcou 
and then to 

(58) k -- 



,-i 



D-a 

We end up with the scaling equation 

(59) — — (aS + xS') + S(DS + xS') + ^- D {DS + x5" - 1) = 0. 

From Eq. (|59jl . we find that the large x asymptotics of S is S(x) ~ x~ a . In a short finite 
time after t co u, it is clear that the large distance behavior of the density profile (r ^> Tq) cannot 
dramatically change. We deduce that the decay of S should match the behavior for time slightly 

2D 

less than t co u for which S(x) ~ x~^+?. Hence the value of a should remain unchanged before 
and after t co u. Finally, we obtain the following exact behaviors for short time after t co u: 

(60) p (t) = ^-{t-t coll )-\ 



D+2 

2 \ 2D D+2 

(61) r Q (t) = (-) (t-t coll ) — .. 



2 

D 



D 

2 



(62) N (t) = fil-) (t-tcou 

We note the remarkable result that the central residual density p(0,t) = S^poft) displays a 
universal behavior just after t co u, a result already obtained in [19J. Moreover, we find that 
No(t) has the same form as the mass found within a sphere of radius ro(t) below t co u, given in 
Eq. flUl). 

Moreover, the scaling function S satisfies 



(63) 



D 2 



(yjj^S + xS^j + S(DS + xS') + px-°(DS + xS' - 1) = 0. 



The constant p is determined by imposing that the large r behavior of s(r, t) (or p(r, t)) exactly 
matches (not simply proportional) that obtained below t co u, which depends on the shape of 
the initial condition as shown in |19j . Equation can be solved implicitly by looking for 
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solutions of the form x D = z[S(x)]. After cumbersome but straightforward calculations, we 
obtain the implicit form 



(64) 



x 



D 



1 + —S(x) 
p 



1 + 



X 



D 



S(x) + 



D 2 



D 

D+2 



which coincides with the implicit solution given in jTHj. Note that S(x) is a function of x D . We 
check that the above result indeed leads to 5(0) = D' 1 , and to the large x asymptotics 



(65) 



S{x) 



D 

D+2 



2D 
X D +' 2 . 



Note finally that for T = 0, N saturates to 1 in a finite time, corresponding to the deter- 
ministic collapse of the outer mass shell initially at r = 1. Indeed, using Gauss' theorem, the 
position of a particle initially at r(t — 0) = 1 satisfies 



(66) 



dr 
~dt 



-(£>-!) 



The position of the outer shell is then 
(67) r(t) 
which vanishes for t en d = -D -1 . 



Dt) 



l/D 



5 Post-collapse dynamics at T > 
5.1 Scaling regime 

In the more general case T ^ 0, we will proceed in a very similar way as in the previous section. 
We define again, 

(68) .(r.O- ^y^ . 
where N still satisfies 

(69) ^ = p N . 
We now obtain 

ds rT1 (d 2 s D + lds\ ( ds ^\ N f ds „ 

(70) — = T — + — +(r— + Ds)s + ^(r— + Ds-po 



dt \dr 2 r dr J \ dr j r D \ dr 

By definition, we have again s(0, t) = p (t)/D. 
We look for a scaling solution of the form 



r 



(71) s(r,t) = p (t)S 

V r ol r J, 

with 5(0) = D^ 1 . As before, we define the King's radius by 

/T \ 1/2 

(72) 

VPo 
12 



For t < t co ii, we had s(r,t) ~ 4Tr~ 2 (or S(x) ~ 4a; -2 ). In a very short time after t co u, this 
property should be preserved, which implies that the post-collapse scaling function should also 
behave as 

(73) S(x) ~ 4x" 2 , 

for large x. Inserting the scaling ansatz in Eq. (|7()jl . we obtain 

( 74 ) A^T? (25 + ar£") = 5" + ^±±S' + SpS + x<S") + ^%^-(DS + xS' - 1). 
2po at x pqTq x u 



Again, this equation should be time independent for scaling to hold, which implies that there 
exist two constants p and k such that 

(75) N = pp r°, 
and 

(tr\ 1 d P° 

(76) = 

with k > 0, as the central residual density is again expected to decrease. Equation (|75jl implies 
that po ~ (t — tcou)' 1 , and then that AT ~ (£ — ico«) D ^ 2_1 - We thus find a power law behavior 
for N , which in order to be compatible with Eq. (J69j) . leads to the universal behavior 

(77) Po (t)= Hl-ljit-tcou)- 1 , 
and then to 

(78) k = — *— . 
K ' D-2 

We end up with the scaling equation 

(79) — *— (2S + xS') + S" + ^^S' + S(DS + xS') + px-°{DS + xS' - 1) = 0, 
D — 2 a; 

where p has to be chosen so that S(x) satisfies the condition of Eq. (|73*jl. Its value will be 
determined numerically in Sec. |U] for D = 3. Note that for small x, the pre-collapse scaling 
function satisfies S(x) — S(0) ~ x 2 , whereas the post-collapse scaling function behaves as 

(80) S(x) - S(0) ~ x D . 

However, contrary to the T = case, S(x) is not purely a function of x D . 

Finally, we find that the weight of the central peak has a universal behavior for short time 
after t coU 

(81) N (t) = p {^-^ D/2 " T D / 2 (t - t coll )^-\ 

Note that No(t) behaves in a very similar manner to the mass within a sphere of radius tq below 
tcou, shown in Eq. (|3*H)l. In addition, comparing Eq. (|ST|) and Eq. (|52*jl. we can define again a 
post-collapse cross-over time between the T^O and T = regimes 

(82) £* — t co u ~ T 15 / 2 , 
which is similar to the definition of Eq. (f4*4*j) . 
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5.2 Large time limit 

Contrary to the T = case, the complete collapse does not take place in a finite time as thermal 
fluctuations always allow for some particle to escape the central strongly attractive potential. 
In order to illustrate this point, and obtain more analytic insight on this matter, we will place 
ourselves in the extreme situation where almost all the mass has collapsed (No ~ 1), and only 
an infinitesimal amount remains in the residual profile. 

In this limit, the residual density p(r,t) satisfies the Fokker-Planck equation 

(83) a ± = T (^l + B^l¥\ + 1 * 



dt \dr 2 r dr J r D 1 dr' 
with boundary condition 

(84) T^ (M ) +p(M)=0 . 

The problem indeed reduces to the study of a very light gas (i.e. with negligible self-interaction) 
of Brownian particles submitted to the gravitational force F = —(GM/r D ^ 1 )e r of a central unit 
mass. Alternatively, this can also be seen as the probability distribution evolution equation of 
a system of two Brownian particles moving in their mutual gravitation field. 

Equation ()8Hj) can be re-expressed as a Schrodinger equation (in imaginary time), thus 
involving a self-adjoint operator (see Appendix [XJ). The large time behavior is dominated by 
the first eigenstate. Coming back to the notation of Eq. (JH3j) . we find that 

(85) p(r,t) ~ e" A V(r), 
where ip satisfies the eigenequation 

(86) -AtKr) = T (> + ^V) + 
and the same boundary condition as p, i.e. 

(87) iy(i)+v(i) = 0. 

The eigenvalue A will also control the large time behavior of po and Nq as Eq. (}6T?j) and Eq. (jHB^) 
both imply that 

(88) i_jv (t)«^ ~e" A *. 

A 

We did not succeed in solving analytically the above eigenequation, and for a given temper- 
ature, this has to be solved numerically. However, in the limit of very small temperature, we 
can apply techniques reminiscent from semiclassical analysis in quantum mechanics (T h). 
We now assume T very small and define </> such that 

ch(r) 

(89) ip(r)= e -—. 

The function h = <$' satisfies the following non-linear first order differential equation 

(90) T (h! + °^±h) + ^ - V? = AT, 
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with the simple boundary condition 



(91) h(l) = l. 

In the limit T — > 0, the term proportional to T in the left-hand side of Eq. (|90|) can a priori 
be discarded leading to 

2XTr D ~ l 

(92) h(r) 



1 + Vl - AXTr 2 ^- 1 ) 

If 4 AT < 1, the above expression is a valid perturbative solution also at r = 1, but cannot 
satisfies the constraint h(l) = 1. Hence, we conclude that in the limit of small temperature 
4AT > 1, so that the above expression is only valid for r not to close to r = 1. The above 
argument also suggests that AT is of order unity and we write 

(93) AT = i + /i 2 . 

To understand how the boundary condition Eq. (|91j) can be in fact satisfied, one has to come 
back to Eq. ()90|) . which for r = 1, shows that h'(l) ~ A ~ T _1 ^> 1. This implies that the term 
Thl cannot be neglected near r = 1, and that h varies in a noticeable way on a length scale 
from 1 of order T. 

This suggests to define 

(94) z{x) = h(l - Tx) 
which satisfies (at order in T) 

(95) -z' + z - z 2 = - + /i 2 , 
and z(0) = 1. This equation has the unique solution 

(96) z(x) = l + ^- 2 ^ X l 

2 2/i + tan(/ix) 

For large x, this function only has a sensible behavior for /i = 0, which shows that 

(97) lim AT = -, 

and that Eq. (J92|) is in fact valid for 1 — r ^> (4AT — 1) — > 0. To leading order, we find 

(98) h(l - Tx) w zq(x) = \ + 1 



2 2 + x 

Equations (}9~2*|) and (|9*Hj) . show that h(x) goes rapidly from 1 to 1/2 in a small region close 
to r = 1 where h varies on the scale T. One can even compute the next correction to AT by 
including the next term of order T in the equation for z. Writing 

(99) z(x) = z (x) + T 1/3 ^(xT 1/3 ), 
we find 

(100) z[ + -z x + z\ - ^-±u = = -c D , u = T^x. 

u 2 _/ ' 
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This is again an eigenvalue problem which selects a unique constant cp, that we could only 
solve numerically. Still, this leads to the non trivial result 

(10!) A = ^ + ^ + - ( T ^°)- 

We now solve the eigenvalue problem ()86|) (JHTj) in the limit of large temperatures T — > +00 
(see also Appendix El). We again perform the change of variables f!89|) and rewrite Eq. (|9Uj) in 
the form 



/ s ,/ D-l, x 1 / h 

(102) + h = A - — I -jr—r - h 



,2 



T \ r 



Then, we expand the solutions of this equation in terms of the small parameter 1/T <C 1. We 
write h = h + ^hi + 5^/12 + ••• and A = Ao + ^Ai + 7^2X2 + •••• To zeroth order, we have 

(103) h' + ^— U = Ao. 

r 

The solution of this equation is ho = X^r/D. Using the boundary condition h${l) = 1, we 
obtain 

(104) A = D, h = r, 
and h n (l) = for n > 0. To first order, we get 

(105) h[ + = \i-^ L - 1 +h 2 . 

Integrating this first order differential equation and using the boundary condition hi(l) = 0, 
we obtain 

D 1 r 3 

(106) hi = — rr - -r 3 -° + , 

K J 2(D + 2) 2 D + 2' 

with 



(107) Ai 



2(D + 2)' 

Hence, the large temperature behavior of the eigenvalue is 

(108) A = fl+ _^l_i ; + ... (T^+oo). 



This expansion can be easily carried out to higher orders but the coefficients are more and more 

JLi _ 477 1 

10 T 700 T 2 



complicated. Restricting ourselves to D = 3, we get A = 3 + y^i — + 0(T 3 ) 



6 Numerical simulations in the canonical ensemble 

In this section, we illustrate the analytical results obtained in the previous section in the 
case of D = 3. Except when specified otherwise, our simulations have been performed at 
T = 1/5 < T c w 0.397..., for which we have obtained t coU w 0.44408.... 
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Figure 1: We plot No(t) for small times (full line). This is compared to 
iVo(t) Theory x [1 + a(t - t coU ) b ] (dashed line), where iVo(t) Theory is given by Eq. (JSTJ with 
ji = 8.38917147..., and a ps 1.7. and b ~ 0.33 are fitting parameters. Note that the valid- 
ity range of this fit goes well beyond the estimated t* with — t co u ~ T D I 2 ~ 0.09. The bottom 
insert illustrates the exponential decay of 1— No(t) ~ e~ A *. The best fit for A leads to A ~ 5.6362 
to be compared to the eigenvalue computed by means of Eq. (j8Tj) . A = 5.6361253.... Finally, 
the top inset illustrates the sensitivity of iVo(t) to the space discretization, which introduces an 
effective cut-off (a factor 4 between each of the 3 curves). Note the small time scale. Even the 
curve corresponding to the coarsest discretization becomes indistinguishable from the others 
for t > 0.448. 
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Figure 2: In the post-collapse regime, we plot p(r,i)/po(t) as a function of the scaling variable 
x = r/r (t). A good data collapse is obtained for central residual densities in the range 
10 3 ~ 10 6 . This is compared to the numerical scaling function computed from Eq. (|79j) (dashed 
line). The insert shows the comparison between this post-collapse scaling function (dashed 
line) and the scaling function below t co u which has been rescaled to have the same value at 
x — 0, preserving the asymptotics: S(x) = (3 + x 2 /4) _1 (see Eq. (jSHJ); full line). Note that the 
post-collapse scaling function is flatter near x = 0, as S(x) — 1/3 ~ x 3 (in D = 3) above t co u 
instead of S(x) — 1/3 ~ x 2 , below t coU . 
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Figure 3: The main plot represents (AT — |) T~ 2 / 3 as a function of T 1//3 (line and squares), 
which should converge to cd =3 = 2.33810741... for T — > according to Eq. (jlUlj) . We find 
a perfect agreement with this value using a quadratic fit (dotted line). Furthermore, this fit 
shows that the slope at T = is in fact equal to —2 ± 2.10~ 4 , suggesting that the next term 
to the expansion of Eq. (jlOlj) is A = ^ + — 2 + in D = 3. In the insert, we plot A as 
a function of T up to T « T c « 0.4. The small temperature analytical result of Eq. (jlOlj) is 
in very good agreement with the numerical data up to T ~ 0.03, whereas the large T estimate 
A(T) = 3 + jqT~ 1 + ... is only qualitatively correct in this range of physical temperatures T < T c . 
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Figure 4: For T = 0.01 (A = 35.074198..., \i = 0.31739877...), we plot h(r) computed numeri- 
cally from Eq. ()90)) (dashed line), and the theoretical expression of Eq. ()92|) . which is valid for 
1 — r ^> /j 2 ~ T 2 / 3 (full line). We also plot the theoretical expression of zq (dashed-dot line; 
see Eq. ()98|)) and the next order perturbation result (dot line; see Eq. (|9*9"j)). which are valid 
in the region 1 — r <C ji 2 ~ T 2 / 3 . The insert is a blow up of the region close to 1. Note how 
h{r) varies by a quantity of order unity as r varies by a quantity of order T = 0.01. We have 
chosen a not too small value for T in order to be able to visualize the two scale regimes on a 
single figure. Both approximations shown in the insert are getting better as T decreases. 
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In order to perform our simulations, we have used a Runge-Kutta algorithm with adaptive 
step in space and time. We call dr the spatial discretization near r = (which we need to 
take very small as the density profile becomes singular at r = 0). An important numerical 
problem arises in the numerical integration of Eq. (|69|h which is crucial in obtaining non zero 
values for N (t). As this equation is a first order differential equation with initial condition 
^o(O) = 0, any naive integration scheme should lead to a strictly vanishing value for N (t) 
for all time, and any dr. Still, when performing this naive numerical integration, we see that 
crossing t co u generates increasing values for M(dr,t), although keeping M(0,t) = N (t) = 
ultimately makes the numerical integration unstable. In order to bypass this problem, we have 
decided to introduce a numerical scheme where Eq. (jHSj) is replaced by 

(109) ^ = pf iVf . 

Nq and Pq are extracted from a fit of M(r, t) to the functional form (we are in D = 3) 

(110) M(r,t) « Nt(t) + ^^-r 3 + a 5 (t)r 5 + a 6 (t)r 6 , 

in a region of a few dr, excluding of course r = 0. This functional form is fully compatible with 
the expected expansion for M{r,t), both below {a§ = 0) and above (a 5 = 0) t co a. We find that 
this numerical scheme allows us to cross smoothly the singularity at t co u. An effective cut-off is 
introduced which effectively depends on dr, and we have checked that the result presented in 
this section are extremely close to the ones that would be obtained in the ideal limit dr — > 0. 
This is illustrated in Fig. 1, where the smoothing effect of our algorithm is shown to act on 
a very small time region after t co u. Even more surprisingly, we find that for sufficiently large 
times (actually very small compared to any physical time scales), our results are essentially 
independent of dr, even for unreasonably large values of dr. We are thus confident that we 
have successfully crossed the collapse singularity. 

In Fig. 1, we plot No(t) for small time which compares well with the universal form of 
Eq. (|81jl. where p = 8.38917147... has been determined so as to ensure the proper behavior of 
S(x) for large x (see Sec. 15. 1|) . We also illustrate the exponential decay of 1 — N (t), with a 
rate in perfect agreement with the value of A extracted from solving numerically the eigenvalue 
problem of Eq. (JSBJ). Finally, we show the effect of the numerical spatial discretization dr near 
r = 0. Satisfactorily enough, the value of N (t) is sensitive to the choice of dr only for very small 
times after the collapse, and we were able to easily reach small enough dr, in order to faithfully 
reproduce the post-collapse singularity. In Fig. 2, we convincingly illustrate the post-collapse 
scaling, and compare the post-collapse scaling function to that obtained analytically below t co u 
(pre-collapse). In Fig. 3, we confirm the validity of our perturbative expansion for A, in the 
limit of small temperature. We compare the value of cd extracted from directly solving the full 
eigenvalue problem to that obtained from Eq. (jlOOj) . finding a perfect agreement. Finally, in 
Fig. 4, we compare the numerical value obtained for h(r) to the different analytical estimates 
given in the preceding section, for T = 0.01. The two important regions 1 — r <C p 2 ~ T 2 ^ 3 
and 1 - r > p 2 ~ T 2 ' 3 can be clearly identified. 



7 Post-collapse in the microcanonical ensemble 

So far, we have only addressed the post-collapse dynamics in the canonical ensemble. In the 
microcanonical ensemble, the dynamical equation have to be supplemented with the strict 
energy conservation condition (see Eq. (J24j) ) . which fixes the global temperature T(t). For this 
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model, it was shown in (HI E3 I2III that below a certain energy (Antonov energy), the system 
collapses with an apparent scaling associated with a max ~ 2.2 for intermediate times (when the 
temperature still increases in a noticeable way) before entering a scaling regime with a = 2, 
identical to that obtained in the canonical ensemble. In the limit t — > t co u, the temperature 
and the potential energy both seem to converge to a finite value preserving a constant energy. 
Closely before the collapse time, the temperature behaves as T{t coU ) — T(t) ~ (t co u — t) 7 with 
7 ~ 1/2. This section addresses the t > t co a time period. 

Assuming a spherical mass density, and after integration by parts, the potential energy W 
can be rewritten in the form (D > 2) 

(m) w{t) = -lfM^ dr . ' 



2 J r D ~ x 2(D-2)' 

We see immediately that as D — 1 > 1, the occurrence of a finite mass iVo(t) 7^ concentrated 
at r = implies an infinite potential energy, hence an infinite temperature. We thus antici- 
pate that the post-collapse dynamics in the microcanonical ensemble is probably an ill-defined 
problem. In this extreme regime, let us try to consider the possible flaws of this model in order 
to describe a consistent dynamics of a reasonable physical self-gravitating system. First, our 
assumption of uniform temperature is certainly not realistic in a system displaying huge den- 
sity contrast, and some alternative approaches are needed to incorporate a spatially dependent 
temperature. This point is certainly crucial and will be addressed in a future work |H2] . Fur- 
thermore, in this regime, a careful physical analysis predicts that this system of self-gravitating 
individual particles should lead to the formation of binaries which is probably beyond the de- 
scription ability of our essentially mean-field approach. In other words, the system may become 
intrinsically heterogeneous, which probably cannot be captured by our continuous model. Fi- 
nally, we can think of other physical effects (degeneracy effects of quantum or dynamical origin, 
finite particle size effects,...) preventing the system from reaching arbitrarily large densities. 
One way to describe such effects consists in introducing a spatial cut-off /iora density cut-off 
of order h~ D . In such a system, the dynamics first follows the pre-collapse dynamics until the 
maximum density is approached. Then, the system will ultimately reach a maximum entropy 
state that we propose to characterize in a simple manner, as in pffi 123]. 

We propose to describe the final state as a "core-halo" structure, which for simplicity we 
modelize as a core of radius h <C 1 and constant density 

(112) Pcore = ^ 

which mimics a regularized central Dirac peak containing a mass Nq. In the region h < r < 1 
stands the halo of constant density 

(113) Phalo - SD(1 _ hD y 

containing the rest of the mass. As h is small, we can compute the potential (or total) energy 
and the entropy only including the relevant leading terms. We find 



(114) E = —T 



£> D 



2 D 2 -4 



m d-i„ d 



+ 1 + ^^ N o-^K 



h D - 2 2 u 2 



0{h 2 



where the first term is the kinetic energy, whereas the entropy (up to irrelevant constants) reads 
(115) S = j InT - N In - (1 - N ) ln(l - N ) + 0(h D ). 
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Figure 5: For E = —0.45 < E c « —0.335, we plot No(t), for different values of h decreasing by 
a factor 2 for each curve from the top one to the bottom one. It is clear that N (t) decreases 
as h decreases. The insert shows the corresponding temperature plots. The temperature T(t) 
increases as h decreases. Note that in the pre-collapse regime, the temperature essentially does 
not depend on h. For the situation considered, it starts from T(0) = 0.1 and culminates at 
T(t coll ) « 0.5. 



For a given small value of h, S has a local maximum at N Q = provided that E > E c (h), 
with lim/^o E c {h) = — £) ^ 4 . Below E c , the sole entropy maximum resides at No satisfying the 
following implicit equation (again in the limit of small h) 

(116) N D ' 



In©) InV 

where the given asymptotics is quantitatively correct only for extremely small values of h. 
Hence, we find that the mass included in the core slowly decreases with the core size h [21], 
resulting in an effective singularity p(r) = — r^8 D (r). Meanwhile, the temperature diverges as 

2 2 1 

(117) T--—W 



D D 2 -Ah D ~ 2 In 2 h 

and leads rapidly and efficiently to a uniform halo. 

In order to relate this result to our actual system, we perform microcanonical post-collapse 
simulations using the same regularization scheme as in the canonical case in order to describe 
the evolution of No(t). In addition to this, we also need to regularize the potential energy which 
is strictly infinite (as well as T) when N ^ 0. Consistently with the previous discussion, we 
introduce a numerical cut-off h by defining 



2 A 2(D-2)' 

In Fig. 5, we plot No(t) and T(t) for different values of the cut-off h. Contrary to the canonical 
case, the post-collapse dynamics is strongly h dependent. We see that in conformity with 
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our result of Eq. (J116)) . the central mass N Q clearly decreases as h —>■ 0. Therefore, in the 
microcanonical case, the physical picture is that when the collapse time t co u is reached, the 
temperature increases rapidly, which leads to the rapid homogenization of the system except 
for a dense and small core, whose mass iV ~ In" 1 T is a decreasing function of the maximum 
temperature reached. This central structure with weak mass and huge binding energy is similar 
to a "binary star" structure in stellar dynamics. Binary formation is the physical process 
that arrests core collapse in globular clusters |34j . This is also the end point of our simple 
microcanonical Brownian model. 

8 Conclusion 

In this paper, we have investigated the post-collapse dynamics of a gas of self-gravitating 
Brownian particles in canonical and microcanonical ensembles. Our results also apply to the 
chemotactic aggregation of bacterial populations in biology. At the collapse time t coU) the 
system develops a singular density profile scaling as p ~ r~ 2 . However, the "central singularity 
contains no mass" , the temperature does not diverge, and the entropy and free energy are finite 
[SI EH 120] • Since this profile is not a maximum entropy (resp. minimum free energy) state 
nor a stationary solution of the Smoluchowski-Poisson system, the collapse continues after t co u. 
This solves the apparent paradox reported in [2IJ. 

In the canonical ensemble, mass accretes progressively at the center of the system and a 
Dirac peak forms by swallowing the surrounding particles. Eventually, the Dirac peak contains 
all the mass. This structure has an infinite free energy F = E — TS —>■ — oo simply because its 
binding energy is infinite. This is therefore the most probable structure in canonical ensemble 
[01 EH 121] ■ In the microcanonical ensemble, the maximum entropy state (at fixed mass and 
energy) consists of a single binary embedded in a hot halo [221 121 CO]- This is precisely what we 
see in our numerical simulations. The temperature increases dramatically above t co u (resulting 
in an almost uniform halo) although the mass contained in the core is weak (but finite). We 
note the "spectacular" fact that almost all the gravitational energy resides in a binary-like 
core with negligible mass. A similar phenomenon is observed in stellar dynamics for globular 
clusters having experienced core collapse PQ. This shows that the microcanonical Smoluchowski- 
Poisson system shares some common properties with kinetic equations usually considered in 
stellar dynamics (Landau- Fokker-Planck equations), despite its greater simplicity. Clearly, a 
major drawback of our microcanonical model is to assume that the temperature uniformizes 
instantaneously, implying an infinite thermal conductivity. We shall relax this simplification in 
a future work [32] • However, the present study is one of the first dynamical study showing the 
formation of Dirac peaks and binary-like structures in systems with gravitational interaction. 

Our Brownian model is based on the existence of a gaseous medium that generates a friction 
force. This situation exists in certain astrophysical models such as the transport of dust particles 
in the solar nebula [12] . Dust particles are submitted to Stokes or Epstein drag. It is clear that 
when the concentration of particles is important (prior to planetesimal formation), self-gravity 
has to be taken into account. Thus, our system of self-gravitating Brownian particles could 
be connected to this astrophysical situation. We just mention this as a possible astrophysical 
application because it is not our present main motivation to make a precise model of dust-gas- 
gravity coupling in protoplanetary disks. However, this problem could be considered in future 
works. 
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A Brownian particle around a "black hole" 

We consider the Brownian motion of a particle subject to the gravitational force —GMr/r 3 
created by a central mass M ("black hole"). We assume that when the particle comes at r = 0, 
it is captured by the central mass. We denote by W(r,t) the density probability of finding the 
particle in r at time t. It is solution of the Fokker-Planck equation 

dW ( r 

(119) -*- = V [TVW + W— 

at \ r 6 

where we have set G = M — R = £ = 1. Let W(r,t) denote a spherically symmetric solution 
of Eq. (|119|) satisfying the boundary conditions 

dW 

(120) T— (l,t)+W(M) = 0, 



(121) ^(r,0) = ^ 



We call 



(122) 3 = -[TVW + W^r 

\ r 6 



the current of probability, i.e. JdSn gives the probability that the particle crosses an element 
of surface dS between t and t + dt (n is a unit vector normal to the element of surface under 
consideration). 

We introduce the probability p(r , t)dt that a particle located initially between r and r +dro 
arrives for the first time at r = between t and t + dt. We have 

f of dW W\ 

(123) p(r , t) = - J ■ dS = Aire 2 T— — + — = AttW(0, t), 



tic 



dr r 



where R € is a ball of radius e —>■ 0. The total probability that the particle initially between tq 
and ro + dro has reached the center of the system between and t is thus 

(124) Q(r ,t)= [ p(r ,t')dt'. 

Jo 

Finally, we average Q(ro,t) over an appropriate range of initial positions in order to get the 
expectation Q(t) that the particle has been captured at time t. 
With the change of variables 

(125) W = ipe^, 

we can transform the Fokker-Planck equation ()119j) into a Schrodinger equation (in imaginary 
time) of the form 

< 126 > w= rA ^-i^ 

A separation of the variables can be effected by the substitution 
(127) ip = (f)(r)e- xt . 
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This transformation reduces the Schrodinger equation to a second order ordinary differential 
equation 



(128) 

with the boundary condition 
(129) 



r 



A 



T 4T 2 r 4 



f 



0. 



ni) + ^0(i) = o. 



We note A n the eigenvalues and (fi n the corresponding eigenf unctions. Since the Schrodinger 
operator H = A — l/4T 2 r 4 is Hermitian, the eigenfunctions form a complete set of orthogonal 
functions for the scalar product 



(fg) = / f(r)g(r)Anr 2 dr. 



(130) 



The system can be furthermore normalized, i.e. {4> n 4> m ) = S nm . Any function f(r) satisfying 
the boundary condition ()129|) can be expanded on this basis, as 



(131) 

In particular, 
(132) 



f(r) = ]T(/0n>0n. 



5(r-r ) 

47rr 2 =2^0n(r o )0 n 



The general solution of the problem f)l 19|) (|120|) can be expressed in the form 
(133) W(r, t) = J2 A„e- A "*e^0 n (r), 

n 

where the coefficients A n are determined by the initial conditions ()121|) . using the expansion 
for the 5-function. We get 



~An,t , 



(134) W(r,t) = e wi ~ ^J^e 

n 

From this expression, we obtain 
(135) 

Then, according to Eq. ()124|) . we have 

(136) Q(r , t) = 4vre"^ ^ ^ r °) lin l 



p(r ,t) = A-ne 2 ^ r oS^ e Xnt (p n (r )\im 



4> n (r)e 7 < 



1 - e~ An * 
1 

A n 



Finally, averaging over the initial conditions, the probability that the particle has been captured 
by the central mass at time t can be expressed as 



(137) 
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where 



(138) 
and 

(139) 



Q n (t) = B n (l - e~ Xnt ) 



1 — ZL 

B n = 4ir—e 2Tr < 



o0 n (r o ) lim 



(r)e 2Tr 



This formally solves the problem. If we consider the large time limit, we just need to determine 
the first eigenvalue Aq(T) of the quantum problem. This has been done analytically in Sec. 15.21 
in the limits T — > and T — > +oo. 

Below, we consider again the high temperature regime where thermal fluctuations prevail 
over gravity but we do not restrict ourselves to the first eigenvalue. To leading order in the 
limit T — > +oo, the Fokker-Planck equation ()119|) reduces to the pure diffusion equation 



(140) 



dW 



,dW 
dr 



However, for consistency (see Sec. I5.2J) . it is necessary to keep the term of order 1/T (arising 
from the gravitational force) in the boundary condition. Hence, we take 



(141) 



dW 1 

— (i,t) + -mi,t) = o. 



The general solution of the diffusion equation ()140|) with the boundary conditions (j!41|) and 
f!121|) can be expressed as 



(142) 

where <fi is solution of 
(143) 



W(r,t) = Y,e- Xnt <p n ( 



r 



V + V + ^0 = 0, 
r 1 



(144) 
Setting 

(145) 



0'(1) + -0(1) = 0. 



x/r, Eqs. (I143j) and (|144j) become 



x" + = o, 



(146) X'(l)= (l-^)x(l). 

Equation (|145|) is readily solved. The eigenvalues can be written \ n = Tx^(T), where x n (T) 
are the solutions of the implicit equation 

(147) tan(x n ) 



1 - -i' 

2T 
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The eigenfunctions are 



(148) Mr)=A n K 

r 

The general solution of the diffusion equation can thus be written 



(149) W(r, t) = e~ Tx ^Al sin(x n r) sin(x n r ), 

with 



^0 n=0 



(150) Al = — ^ 

27r{x n — smx n cosx n ) 

If we consider the pure diffusion of a particle in a box, the boundary condition ()144|) reduces 
to <p'(l) = and the solutions of the implicit equation 

(151) tan(x n ) = x n . 

In particular, xq = 0. This implies that the probability W(r,t) converges for large times to 
a uniform profile W(r, +00) = 3/Att which is indeed solution of the diffusion equation in a 
box. If gravity is taken into account, its first order effect (in the limit T — > +00) is to change 
the boundary condition to Eq. ()144j) . It is as if we had a diffusion across the box [23 IT2] 
although the true physical process is a capture by the central mass. The eigenvalues are now 
determined by Eq. ()147j) . The hardly modified (to first order) with respect to the 

preceding problem but xq is now different from zero. To first order, we find that Xq = 3/T so 
that Ao = 3 in agreement with the result of Sec. 15.21 We also note that Aq = T /Air while A n>0 
are independent on T (to leading order) and given by Eqs. (|15U|) and (|151|) . 

Using Eqs. ()123|) and ()124|) . the probability that the particle has been captured by the 
central mass at time t is given by 

(152) Qffl^^lz^^-M 



1 n=0 Xn r ° 



If we average over initial conditions with the weight 3r$ (uniform distribution), we find to 
leading order in T _1 that 



(153) = 0, for n > 0, 

r 



(154) Sin(Xoro) = x . 

Hence, the modes n > cancel out. Therefore, in the high temperature regime, the probability 
of capture is given by 

(155) Q{t) = 1 - e~ 3t , 
for all times. 
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The case D = 2 can be treated by a similar method. Instead of Eqs. (|149j) (|15(Jj) and (J147)) . 

we get 

+00 

(156) W(r,t) = J2 e ~ TxltA2 n J oMJ (x n r ), 

n=0 

(157) A\ 



(158) 



n[Jl{x n ) + J 2 (x n )]' 

XnJli^Xn) 1 

Jo(x n ) T" 



where J n is the Bessel function of order n. For the pure diffusion process, x n = a\ n are the 
zeros of J\. If gravity is taken into account, then x n> \ ~ a± n while Xq = 2/T establishing 
Ao = 2. The probability that the particle has been captured by the central mass at time t is 
given by 

Ott + °° 1 _ P -Tx 2 n t 

(159) Q(t) = — % AlJ (x n r ). 

1 n=0 %n 

If we average over the initial conditions with a weight 2r$ (uniform distribution), we get 
Jo(x n r ) = if n > and Jo(x r ) = 1. Therefore, in the high temperature regime, the 
probability of capture is given, for all times, by 

(160) Q(t) = 1 - e~ 2t . 
Finally, for D — 1, we obtain 

+00 

(161) W(r, t) = ^] e~ Tx2nt A 2 n cos(x n r) cos(a; n r ), 

n=0 



(162) Al 



2 ^ 

sin(2a: n )i 



[1 



2,i' n 



(163) x„tan(x„) = ^. 

For the pure diffusion process, x n = nn. If gravity is taken into account, then x n> \ ~ nil while 
Xq = 1/T establishing A = 1. The probability that the particle has been captured by the 
central mass at time t is given by 

(164) Q(t) = -J2 2 K^s(x n r ). 

1 n=0 Xn 



If we average over the initial conditions with a weight 1 (uniform distribution), we get cos(x n r ) = 
if n > and cos(xor ) = 1. Therefore, in the high temperature regime, the probability of 
capture is given, for all times, by 

(165) Q(t) = 1 - e~*. 
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